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Abstract — This paper presents a boundary-integral equation 
(BIE) method for the calculation of poloidal axisymmetric mag- 
netic fields applicable in a wide range of ac frequencies. The 
method is based on the vector potential formulation and it 
uses the Green's functions of Laplace and Helmholtz equations 
for the exterior and interior of conductors, respectively. The 
work is particularly focused on a calculation of axisymmetric 
Green's function for the Helmholtz equation which is both 
simpler and more accurate compared to previous approaches. 
Three different approaches are used for calculation of the Green's 
function depending on the parameter range. For low and high 
dimensionless ac frequencies we use a power series expansion 
in terms of elliptical integrals and an asymptotic series in 
terms of modified Bessel functions of second kind, respectively. 
For the intermediate frequency range, Gauss-Chebyshev-Lobatto 
quadratures are used. The method is verified by comparing with 
the analytical solution for a sphere in a uniform external ac field. 
The application of the method is demonstrated for a composite 
model inductor containing an external secondary circuit. 

Index Terms — Integral equations. Green function, Helmholtz 
equations. Boundary-element method. Electrical engineering 
computing 



I. Introduction 

CALCULATION of alternating magnetic fields and the 
associated eddy currents is important for the design of 
various electrical machines and the magnetic field inductors 
used for heating, melting, stirring, shaping or levitation of 
metallic or semiconducting materials. Although the distri- 
bution of electromagnetic fields is, in principle, completely 
described by the Maxwell equations, only in very few simple 
cases these equations can be solved analytically. Usually a 
numerical approach is needed. 

Most of the approaches used for the solution of electromag- 
netic field problems are based on finite difference (FDM) or 
finite element methods (FEM). The main advantage of these 
methods is their capability to deal with complex geometrical 
configurations usually encountered in practical applications. 
However, these methods involve a solution for the fields in 
the free space which is often unbounded and, thus, may 
require considerable additional computer resources as well as 
a special numerical treatment at the outer open boundary [1]. 
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On the other hand, in many applications such as, for example, 
electromagnetic heating and hardening of workpieces or the 
stirring of molten metals, only the eddy currents and magnetic 
fields in the conducting medium are needed. There are several 
approaches avoiding the solution of the magnetic field in the 
free space. A first kind of those approaches uses the Biot- 
Savart law, to reduce the problem to a volume integral equation 
[2], [3], [4]. A second type of approaches combine FEM for 
the solution of the corresponding partial differential equation 
(PDF) inside the conductor with a boundary element method 
(BEM) based on the second Green's theorem to represent 
the field in the free space as an integral of the field values 
and its gradient over the surface of the conductor [5], [6]. In 
case of a thin skin-layer this approach reduces to a boundary 
integral equation [7]. On the other hand, there are approaches 
using the boundary impedance condition to approximate the 
field distribution in the conductor at small skin-depth in 
combination with a FEM solution in the free space [8]. 

Our approach is to reduce the problem to a boundary 
integral equation not only for the free space but also for the 
interior of the conductor that would be applicable regardless 
of the relative thickness of the skin-layer This is possible 
for a conductor with uniform electric and magnetic properties 
when the field distribution inside the conductor is described by 
linear PDEs with constant coefficients admitting an analytic 
fundamental solution, i.e., the Green's function. Advantage 
of this approach is a consideration of the conductor surface 
only. Thus, the dimensionality of the problem is reduced 
by one that renders this method particularly suited for the 
analysis of complicated geometrical configurations. On the 
other hand, this geometrical simplification comes at the price 
of an increased algebraic and numeric complexity due to 
the calculation of the axisymmetric Green's function. Similar 
approaches have already been considered for 2D [9], [10] and 
3D [11] problems which both are considerably algebraically 
simpler than the axisymmetric problem considered here. There 
is an analytic solution for the Green's function for the 2D case 
and a point-source Green's function is used in the 3D case 
while there is no simple analytic solution for Green's function 
in the axisymmetric case. The axisymmetric case has been 
addressed by Fawzi et al. [19] who derived boundary integral 
equations for the transverse magnetic (TM) mode in terms 
of the azimuthal electrical and tangential magnetic fields in 
the full electrodynamic formulation including the displacement 
current. The same problem has been revisited in Refs. [20], 
[21] in a quasi-static approximation. Our approach differs 
from the previous ones by a more exact calculation of the 
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Green's function using a combination of analytic, asymptotic 
and numeric methods. 

The paper is organized as follows. In Section II, problem 
formulation and basic equations are given. The boundary- 
integral equation derivation and the calculation of the Green's 
function for the azimuthal component of the vector potential 
is presented in Section III. In Section IV, we describe the 
numerical implementation of the method and give several 
application examples. Finally, summary and conclusions are 
given in Section V. 

II. Problem formulation and basic equations 

Consider an axisymmetric body of a characteristic size Rq 
at rest having a uniform electrical conductivity a placed in 
an axisymmetric external ac magnetic field with induction B 
alternating harmonically in time with the angular frequency 
ut. Searching for magnetic and electric fields in terms of 
vector and scalar potentials as A (r, t) ~ 3t [A{r)e^^*] and 
<i> (r,t) = [$(r)e'"*] , where A(r) and $(r) are generally 
complex axisymmetric amplitudes, leads to the governing 
equation: 



ILUA - 



-V X V X A = -V$. 



In Eq. Q the gradient of the scalar potential V$ plays the 
role of a source term with respect to the vector potential. As it 
will be shown later, this source term can be used to specify an 
externally applied ac voltage to an axisymmetric coil system. 
For the following it is advantageous to use the transformation 



A' 



(2) 



that allows us to remove the source term from Eq. Q by 
including it into the vector potential. Then the equation for 
A' satisfying the Coulomb gauge V • A' = can be written 
as 

V^A' - X^A' = 0, (3) 

where = ioj with uj = noaujRQ being the dimensionless 
frequency. Henceforth all quantities and differential operators 
are supposed to be nondimensionalized by using the corre- 
sponding characteristic length and vector potential scales, i?o 
and Aq , where the latter will be specified in the following for 
each particular problem. The advantage of Eq. (|3} compared 
to its nonuniform counterpart (Eq.Q is that the solution of the 
former can straightforwardly be written as a surface integral 
which is the aim of the next section. 

III. Boundary-integral equation for the vector 

POTENTIAL 

Using second Green's vector theorem the solution of Eq. 
(|3} satisfying the Coulomb gauge can be written as: 

A'(r) = ^ / [VG^(r' - r) x n' x A'(r') 

^ s 

- VG^(r'-r) (n'- A'(r')) 

- G^{r' - r)n' x V x A'{r')] Sr' (4) 



where G'^{r) 



cxp(-A|r|) 
\r\ 



from that of Huang et al. [20] who uses a scalar counterpart of 
the second Green identity which is not correct in general but 
leads to the right result in the special case of an axisymmetric 
and purely azimuthal vector potential which satisfies the 
Coulomb gauge straightforwardly. To find A' at the point r 
inside the volume enclosed by surface S we need the values 
of both a! and n' x V x A' on S. Usually both these values 
are unknown and two vector equations are needed to find 
them. The first of the equations is obtained by approaching 
the observation point r to the boundary 5* that results in: 

/ [VG^(r' - r) X n' X A'(r') - G^(r' - r) 
Js 

n' X V X A'(r')] d^r' - 27rc(r)A'(r) ^ 0, (5) 

where c(r) is a geometrical parameter which is equal to unity 
for a smooth surface [12]. The second equation is obtained by 
considering the nonconducting space outside the body where 
the distribution of the vector potential is governed by a Laplace 
equation. The corresponding equation takes the form 

/ [VG°(r' -r)xnx A-G°{r' -r)nxV X A] (fr' 
Js 

+ 27rc(r)A(r) = 0, (6) 

where the sign difference at the second term is because of 
n is directed inwards with respect to the region outside the 
body. Equation (|5} can now be represented back in terms of 
the original vector potential A and the imposed gradient of 
the scalar potential by inverting the transformation given by 
Eq. (|3: A' = A-iuj-^V<P. 

In the following we focus on the case of a purely az- 
imuthal and axisymmetric vector potential A{r) = e^A{r, z) 
depending only on the radius r and the axial coordinate z in 
a cylindrical system of coordinates. For the gradient of the 
scalar potential V<I> to be purely azimuthal and axisymmetric, 
<I> can be a function of the azimuthal angle ip only: $ = ^{(fi)- 
Then V$ = e^i|| with ff = *o = const because of 
the axisymmetry. Further note that for axisymmetric bodies 
with simply connected shapes including the symmetry axis, 
$0 must be zero for V<1> to be limited at the symmetry axis 
r = 0. However, $o may be nonzero for toroidal bodies which 
are not intersected by the symmetry axis. For such bodies, like 
coils, $0 may be used to specify the externally applied voltage 
U driving the current as §^ = 5^- Alternatively, $0 may be 
determined in the course of solution when the total current 
rather than the voltage is specified on the coil. Note that our 
treatment of the source term is more mathematically rigorous 
compared to Ref. [20]. 

Substituting such a purely azimuthal vector potential into 
Eq. (|6j and performing the integration along the azimuthal 
angle ip we obtain after some transformations an equation 
defining A outside the conducting body 

1 



Air) 



in 



djr'Air')) 

dn' ^^^"^'"^^ 

dn' 



(7) 



is the Green's function of the 



scalar Heknholtz equation. Note that our approach here differs 



where the integral is now evaluated along the contour L 
forming the conducting body of rotation. The Green's function 
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for the azimuthal component of the vector potential entering 
the above equation is 



Jo 

2 sin^ - 1 



2k 



Tr/2 



r'r Jo y/l- sin^ (p 



dip 



Ak 



K{k) - E{k) K{k) 



which is the vector potential of a circular current loop divided 
by r' [13] presented in terms of the complete elliptical inte- 
grals of the first and second kind, K{k) and E{k), respectively, 
of the modulus k = 2^/ , , , yf.^, , — ^ [141. Thus, the 
Green's function like its gradient for the azimuthal component 
of the Laplace equation is obtained analytically. 

The azimuthal component of the vector potential inside the 
conducting body can be obtained in a similar way as outside 
by using the corresponding Green's function with A 7^ 



where 



2k 



dn' ^^^"^'"^^ 
dn' 



r e'^G^r' ~ r)dp' 
Jo 



2 sin p - 1 



dlr'l 



(9) 



7r/2 



^1 — fc^ sin^ p 



rr Ja 



X exp y—ny 1 — fc2 sin pj dp, (10) 

and n — 2A\/rV/fc. In contrast to the previous case with 
A = 0, the last integral cannot be evaluated analytically. For 
|k| <C 1, corresponding to low frequencies, the exponential 
function in dlOt may be expanded into a power series of k : 



2k ^ (-k)" 



ri=0 



1 



(11) 



where /„ 



n ■ 
n 



21 
21 



1 



Jo 



7r/2 



I 



(1 - fc^ sin^ ip) ^ dp 
0,1,2,.., and /; = 

?or odd n the theory of elliptical integrals [15] yields the 
following recursion 



^'+2 - 21 



with Iq 



1 + 2 

f and If ^ f (2 - A:2). Derivative of this 
recursion with respect to fc^ leads to a similar recursion for 
1° . Similarly, for even indices one obtains: 

2^ + 2(2-fc2)/r^,-|±l(i-fcVr 



1+2 



21 



21 



with /q = K{k) and Ii — E{k). Series (II 1> is summed until 



< 10 that ensures a relative error less than 10 for 



At high frequencies, when |k| ^ 1, jlO> is dominated by the 
maximum of the exponential function about the point </? = f 
and it is possible to evaluate it asymptotically by the Laplace 
method [16]. Substitution of cosp ~ t in Eq. dlOt results in 



2(3 f exp (-sv/TTW) (1 - 2^) 



dt. 



where s = k\J\ — k^ and /3 = 



r'r Jq Jl + (3H^ ' 

(12) 

Since the dominating 
contribution in the above integral results from the vicinity of 
t = we can expand = y,°" n ^111+1/^^2™ ^ j^jf^ 

the upper limit of integration to infinity 



G^(r,r') 



r'r Jo 



exp (— s cosh x) 1 — 2 




m=0 
00 



2 m 



Er(rn + 1/2) / sinhx\ 
— I — — I dx 



= E 



r(TO + 1/2) 



'r'r — ^ 0rrn!/32 

m— 



2 

/?2 



where we have made the additional substitution t = 



Wi), (13) 

sinh X 



The integrals in the above relation 

Im — exp (— s coshx) sinh^™ xdx 
Jo 

r(m+ 1/2) /2' 



Km{s), 



defined in terms of the modified Bessel function of the second 
kind of order m, Km{s), [14], can efficiently be calculated for 
m > 1 by the following recursion: 



= (2m + 1) (2m/™ + (27n - l)/™_i) /s^. 



The gradient of can be found in a similar way by using 



the relation 



dl„ 

ds 



im+T which follows from the properties 



of Bessel functions [14]. 

There is an additional range of parameters where the power 
series solution given by Eq. Jl 1> is not applicable because 
|k| is large, while the asymptotic approximation ( I13l l does not 
work because k is small and the exponential function under 
the integral (I10> varies weakly along the angle p without 
having a pronounced maximum. In this case, one could expand 
the sub-integral function in Eq. (I10> in a power series of 
k^. As easy to see, this would result in the power series of 
sin^ (p which can in principle be integrated analytically term 
by term. On the other hand, such polynomials can efficiently 
be integrated by Gauss-Chebyshev-Lobatto quadratures. Thus, 
instead of expanding the integral (I10> in a power series of 
small fc2 and then integrating analytically term by term, we 
apply a Gauss-Chebyshev-Lobatto quadrature [14] directly to 
the integral (I12> . 

To summarize, three different approaches are used for the 
evaluation of the Green's function and its gradient for A 7^ 
within the following parameter ranges defined in terms of fc^ 
and \k\ which actually specify the integral in ( I10> . First, for 
sufficiently small \k\ < 5k^ we use the power series expansion 
(II 1> . Second, for the intermediate range Sfc^ < |k| < 35k^^ a 
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Gauss-Chebyshev-Lobatto quadrature with M — 30 + 120fc^ 
number of points is used where the number of points is 
increased as A: 1 in order to ensure sufficient accuracy in the 
vicinity of the singularity at fc 1. In addition, for > 0.98 
we subtract the singularity as the zero-frequency Green's 
function which can be integrated analytically whereas the rest 
is integrated numerically as described above. For |k| > 35k^^ 
the first five terms of the asymptotic series ( I13> are used. 
The ranges of applicability of different approximations and 
the number of the quadrature points are found numerically 
and they ensure the relative error of the Green's function and 
its gradient to be below 10"'' for k < 0.999. 

Two coupled boundary-integral equations are obtained from 
Eqs. and (|9} by taking the observation point r to the surface 
contour L : 



r' dn' 

- 27rc(r)*(r) = 0; (14) 



dn' 



^'{r')d{r'rG^^{r,r')) 



r' dn' 

+ 27rc(r)«''(r) 0, (15) 



d\r'\ 



where the unknown functions to be found along L are 



^{r) = rA{r) and 



(9*(r) 

dn 



while ^'{r) = ^{r) - iw^^^o 



for the interior involves an additional constant $o defining 
the azimuthal gradient of the electrostatic potential which, 
as discussed above, may be non-zero for the conductors not 
intersected by the symmetry axis. For such conductors the 
geometrical parameter in Eqs. (1141 115> can be determined as 



ld{r'rG%{ry)) 



d\r'\ 



(16) 



which follows from the requirement for Eq. M5\ to be satisfied 
by a constant in the limit of A ^ similarly to its PDF 
counterpart (|3}. 

IV. Numerical implementation and examples of 

APPLICATION 



The system of two coupled boundary-integral Fqs. ( I14II15> 
can be solved numerically by the boundary element method 
[12]. For this purpose each line Lk forming a closed surface 
of a part of the conducting body of rotation, which may 
be simply or multiply connected, is approximated by N 
rectilinear segments with endpoints pointed by radius vectors 
.^N + 1. The integrals in Fqs. 



r^, I = 1,..,A^ + 1. The integrals in Fqs. (I14I16> along 
each contour are replaced by the sums over the corresponding 
boundary elements where the integrals over each boundary 
element are approximated by four-point Gauss quadratures 
[17]. When the observation and integration points coincide 
there is a logarithmic singularity in the Green's function which 
is subtracted and integrated analytically over the corresponding 
element. In the simplest case, the unknown functions are 
considered to be constant within each element that results in 
2N unknown quantities which are the constant values of '^{r) 
and 



On 



in each element. Upon evaluation of both Fqs. ( 1141 
I15> at the midpoint of each element we obtain a system of 2N 



complex linear equations. For a typical number of unknowns 
of about several hundreds this problem can straightforwardly 
be solved by an LU decomposition. 

In the following, we consider two simple examples of appli- 
cation of the method. The first example is a conducting sphere 
of radius i?o in a uniform external ac magnetic field with 
induction amplitude Bq. In this case, the contour encircling 
the whole free space in Fq. MAI may be considered to consist 
of two contours where L encloses the sphere while the second 
one encloses some remote inductor creating a uniform field 
with 5*0 (t") = r^/2 that corresponds to the vector potential 
scaled by = BqEq. 

Thus Fq. il4i for the outer surface of the sphere takes the 
form 



a*(r') 
dn' 



^{r')d{r'rG%{r,r')) 



dn' 



d\r' 



+ 27rc(r)5'(r) = -47rv['o(r), 

whereas the corresponding Fq. dlSI l for the inner surface 
remains unchanged. The distributions of ^{r) and ^^q^^ 
calculated with = 30 constant surface elements are seen 
in Fig. [2 to be in good agreement with the corresponding 



analytical solutions [18]: 

= f 1 _ 1 Mil 
|r| = l V 2jo(a;) 



|r|=i 



3 ji(x) 
2 xjQ(x) 



sin [9) and 



fi\T?[9), where x = \Jio/i^ 









dn r 



dn 

the poloidal angle, and jn{x) is the spherical Bessel function 
of order n [14]. 

An additional quantity which can be used for verification 
of the method is the total dissipated power defined in terms 
of dimensionless surface quantities as 



P 



where the asterisk denotes the complex conjugate and the 
power is scaled by Pq = ^"^^^ . Comparison of numerical 
and exact solutions of total power for a sphere in a uniform 
ac magnetic field plotted in Fig. |2ja) shows that 30 constant 
boundary elements ensure a relative error below a few per cent 
for the dimensionless frequency up to 10'^. For comparison 
we show also the relative error of the solution resulting from 
purely numerical calculation of the Green's function and its 
gradient, as in Ref. [20], with 64 and 128 Gauss-Chebyshev- 
Lobatto quadrature points that results in a significantly lower 
accuracy at both low and high frequencies. As seen in Fig. 
|2jb), the accuracy decreases at high frequencies where a 
larger number of boundary elements is required. Note that 
the relatively slow convergence rate of about ~ is due to 
the low accuracy of constant boundary elements used in this 
example. 

As a next example we consider a model inductor consisting 
of two coaxial mirror-symmetric rings of trapezoidal cross- 
section as shown in Fig.|3]related to the crystal growth appli- 
cation by the floating zone technique [24]. The upper, primary, 
ring defined by the contour Lq is connected to a power source 
supplying ac current / = Iocos{ujt). The current in the 
lower, secondary, ring Li, which is short-circuited through 
an additional impedance Z2, is induced only by the magnetic 
field of the upper ring. In this case, we have two additional 
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poloidal angle, 9 (7r) poloidal angle, 6 (It) 




poloidal angle, 9 (Jl) poloidal angle. 6 (7t) 



Fig. 1 . Comparison of numerical (dots) and analytical (solid curves) solutions for a sphere in a uniform ac magnetic field; real (left) and imaginary (right) 
parts of * (top) and §^ (bottom) at the surface of the sphere versus the poloidal angle for various dimensionless ac frequencies. 




10"' 10° lo' 10^ lo' 16 32 64 128 



Dimensionless frequency, 01 Number of elements, N 

(a) (b) 

Fig. 2. Comparison of the numerically obtained dimensionless dissipated power with the exact analytic solution for a sphere in a uniform ac field: power 
and the relative errors resulting from combined analytical/numerical and purely numerical calculations of the Green's function with M = 64 and M = 128 
Gauss-Chebyshev-Lobatto quadrature points versus the dimensionless frequency for iV = 30 boundary elements (a) and relative error versus the number of 
BEM at various dimensionless ac frequencies (b). 
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unknown quantities, the azimuthal gradients of the electrostatic 
potential $o and $i in the primary and secondary inductors. 
Respectively, we have two additional equations for the circuits 
of primary and secondary rings Iq = J^^ ^ ^'olT'' ^ ™d 

Ili r ^ alr^ 1^1 ~ 27r$i/Z2 = 0, which are discretized and 
solved as described above. Our goal here is to choose the 
additional impedance Z2 so that to have the current induced in 
the secondary ring of the same amplitude but delayed in phase 
by 90 degrees with respect to the current in the primary ring. 
The corresponding magnetic flux lines calculated using 150 
equally-sized boundary elements in each ring are plotted in 
Fig-EJfor various dimensionless frequencies. It has to be noted 
that the problem tends to be ill-conditioned at cZ; = and, thus, 
special care is necessary for calculation of the geometrical 
parameter c(r) defined by Eq. (I16> in order to obtain accurate 
solutions at low frequencies ui 1 which, however, are not 
very important for practical applications. The dimensionless 
impedances of primary and secondary circuits are plotted in 
Fig. Ill a) for the secondary current of the same amplitude but 
tt/2 phase lag relative to the primary one. As seen, the results 
of the present approach are in good agreement with those of 
the boundary impedance condition (BIC) approximation which 
becomes applicable at sufficiently high frequencies (w > 10^). 
The results shown in Fig. |4la) imply that for (D > 30 the 
impedance of the secondary ring can be compensated by an 
additional impedance Z2 ~ R + containing active and 
capacitative components denoted by R and C, respectively. 
The relative amplitudes of the secondary current are plotted 
in Fig. 01b) versus the dimensionless frequency for various 
impedances added in the circuit of the secondary coil. The 
solid curve corresponds to the short-circuited secondary ring. 
As seen, the capacitance determines the resonance frequency 
at which the amplitude of the secondary current attains a 
maximum while its phase becomes delayed by tt/2 with 
respect to that of the primary current. The resistance is added 
to balance the amplitude of the secondary current with that of 
the primary one at the resonance frequency. 

V. Summary and conclusions 

We have presented a boundary-integral equation method 
for the calculation of poloidal axisymmetric magnetic fields 
applicable in a wide range of ac frequencies. The method is 
based on the vector potential formulation and uses Green's 
functions for Laplace and Helmholtz equations for the exterior 
and interior of conductors. Particular attention was paid to the 
calculation of Green's function for the Helmholtz equation 
which underlies our approach. In contrast to the Laplace 
equation, there is no simple analytic solution for the ax- 
isymmetric Green's function of the Helmholtz equation. Thus, 
the corresponding function as well as its gradient has to be 
calculated numerically that is done by three different ap- 
proaches depending on the parameter range. For low and high 
dimensionless frequencies we use power series expansions in 
terms of elliptical integrals and asymptotic series in terms of 
modified Bessel functions, respectively. For the intermediate 
frequency range, Gauss-Chebyshev-Lobatto quadratures are 
used. 



Our way of calculation of the Green's function differs 
considerably from previous approaches. Note that, on the one 
hand, our derivation of the axisymmetric Green's function is 
more straightforward and leads to considerably simpler ana- 
lytic expressions compared to the Fourier series representation 
in terms of Bessel/Hankel and Legendre functions obtained 
by variable separation in spherical coordinates [21] or to the 
Fourier/Bessel integrals obtained by the corresponding integral 
transforms in the cylindrical coordinates [22], [23]. Fourier se- 
ries and the corresponding integrals are computationally more 
expensive because they contain products of special functions 
of a varying argument whereas the power and asymptotic 
series in our case contain only a single special function of 
a fixed argument and varying order which can efficiently be 
calculated using recursion. Moreover, the selective calculation 
of the Green's function by either numerical quadratures, power 
or asymptotic series depending on its argument provides the 
best convergence in each parameter range and, thus, it is 
obviously more efficient numerically than a general Fourier 
series or corresponding integrals. On the other hand, the way 
in which we calculate the Green's function differs significantly 
from the approach of Huang et al. [20] who evaluate integral 
(I10> numerically. Although before integration they subtract 
the Green's function for the Laplace equation in order to 
remove the singularity from the integrand, the discontinuities 
still remain in derivatives and deteriorate the accuracy of the 
numerical integration when the observation point approaches 
the contour of integration. Additional difficulties with numer- 
ical integration arise at high dimensionless frequencies when 
the exponential function in the integrand MQ\ decays in a fast 
and oscillatory way. 

We avoid these problems by calculating the Green's function 
analytically in the form of power and asymptotic series for low 
and high frequencies, respectively. 

The method was verified by comparison with the analytic 
solution for a sphere in a uniform ac magnetic field. In 
addition, the performance of the method was demonstrated 
for a composite model inductor supplied with a current of 
fixed amplitude and containing a secondary coil with an 
external circuit. In this case, the results were checked by 
comparison with an approximate solution obtained by the 
boundary impedance condition which becomes applicable at 
sufficiently high ac frequencies. The accuracy of the numerical 
solution deteriorates at very low frequencies where an increase 
of the number of boundary elements is necessary to obtain a 
smooth distribution of the magnetic field component shifted 
in phase by 7r/2 with respect to the applied potential. 

The proposed method is well suited for the numerical 
calculation of axisymmetric poloidal magnetic field inductors 
of complicated geometrical configurations at intermediate ac 
frequencies because it requires only the surface but no spatial 
discretization. 
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